home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_8 / grads.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  2.6 KB  |  109 lines  |  [MATF/MATL]
  1. function [P0,y0,h,err,P] = grads(Fn,Gn,P0,max1,delta,epsilon,show)
  2. % [P0,y0,P] = grads(Fn,Gn,P0,max1,delta,epsilon,show)
  3. % Gradient search for a minimum.
  4. % Fn is the function, input.
  5. % Gn is the gradient for Fn, input.
  6. % P0 is the starting point, input.
  7. % max1 is the maximum number of iterations, input
  8. % delta   is the tolerance, input;
  9. % epsilon is the tolerance, input;
  10. % show  if show==1 the iterations are displayed, input.
  11. % P0 is the point for the minimum, output.
  12. % y0 is the function value  Fn(V0), output.
  13. % dV is the size of the final simplex, output.
  14. % P  is a vector containing the iterations, output.
  15. if nargin==5, show = 0; end
  16. [mm n] = size(P0);
  17. maxj = 20;
  18. big = 1e8;
  19. h = 1;
  20. len = norm(P0);
  21. y0 = Fn(P0);
  22. if (len>1e4), h = len/1e4; end
  23. err = 1;
  24. cnt = 0;
  25. cond = 0;
  26. while (cnt<max1 & cond~=5 & (h>delta | err>epsilon))
  27.   S = feval(Gn,P0); echo off;
  28.   P1 = P0 + h*S;
  29.   P2 = P0 + 2*h*S;
  30.   y1 = feval(Fn,P1);
  31.   y2 = feval(Fn,P2);
  32.   cond = 0;
  33.   j = 0;
  34.   while (j<maxj & cond==0)
  35.     len = norm(P0);
  36.     if (y0<y1),
  37.       P2 = P1;
  38.       y2 = y1;
  39.       h = h/2;
  40.       P1 = P0 + h*S;
  41.       y1 = feval(Fn,P1);
  42.     else
  43.       if (y2<y1),
  44.         P1 = P2;
  45.         y1 = y2;
  46.         h = 2*h;
  47.         P2 = P0 + 2*h*S;
  48.         y2 = feval(Fn,P2);
  49.       else
  50.         cond = -1;
  51.       end
  52.     end
  53.     j = j+1;
  54.     if (h<delta), cond=1; end
  55.     if (abs(h)>big | len>big), cond=5; end
  56.   end
  57.   if (cond==5),
  58.     Pmin = P1;
  59.     ymin = y1;
  60.   else
  61.     d = 4*y1 - 2*y0 - 2*y2;     % Start of a long block:
  62.     if (d<0),
  63.       hmin = h*(4*y1-3*y0-y2)/d;
  64.     else
  65.       cond = 4;
  66.       hmin = h/3;
  67.     end
  68.     Pmin = P0 + hmin*S;
  69.     ymin = feval(Fn,Pmin);
  70.     h0 = abs(hmin);
  71.     h1 = abs(hmin-h);
  72.     h2 = abs(hmin-2*h);
  73.     if (h0<h), h = h0; end
  74.     if (h1<h), h = h1; end
  75.     if (h2<h), h = h2; end
  76.     if (h==0), h = hmin; end
  77.     if (h<delta), cond=1; end
  78.     e0 = abs(y0-ymin);
  79.     e1 = abs(y1-ymin);
  80.     e2 = abs(y2-ymin);
  81.     if (e0~=0 & e0<err), err = e0; end
  82.     if (e1~=0 & e1<err), err = e1; end
  83.     if (e2~=0 & e2<err), err = e2; end
  84.     if (e0==0 & e1==0 & e2==0), err = 0; end
  85.     if (err<epsilon), cond=2; end
  86.     if (cond==2 & h<delta), cond=3; end
  87.   end     % End of the long block.
  88.   cnt = cnt+1;
  89.   P0 = Pmin;
  90.   y0 = ymin;
  91.   home;
  92.   if cnt==1,
  93.     diary output,...
  94.     disp('The results for a Gradient search.'),...
  95.     disp('      p         q       f(p,q)'),...
  96.     diary off,clc; 
  97.   end;
  98.   if show==1,
  99.     Mx1 = 'Gradient search iteration No. ';
  100.     Mx2 = '      p         q       f(p,q)';
  101.     disp([Mx1,int2str(cnt)]),disp(Mx2),...
  102.     diary output,disp([P0,y0]),diary off;
  103.     pause(0.5);
  104.     plot(P0(1),P0(2),'or');
  105.   end
  106. end
  107.  
  108.  
  109.